# Replication code for Dotan Haim, Matthew Nanes, and Michael Davidson, "Family Matters: The Double-Edged Sword of Police-Community Connections." Journal of Politics.
# Code intended for 'R' Statistical software. Corresponding author Dotan Haim, dhaim@fsu.edu

###############################
######### Preparatory #########
###############################

  # Install packages
  install.packages("ggplot2")  
  install.packages("ggpubr")      
  install.packages("data.table")
  install.packages("readstata13")
  install.packages("igraph")
  install.packages("gridExtra")
  install.packages("RColorBrewer")

  # Load packages
  suppressMessages(library("ggplot2"))  
  suppressMessages(library("ggpubr"))      
  suppressMessages(library("data.table"))
  suppressMessages(library("readstata13"))
  suppressMessages(library("igraph"))
  suppressMessages(library("gridExtra"))
  suppressMessages(library("RColorBrewer"))

  # Working Directory
  setwd("FOLDER WHERE YOU UNZIPPED THE REPLICATION FILE")

###########################################
######### Figure 1: Network Plots #########
###########################################

  # ===========================================================
  # Prepare Network Attributes for Plotting
  # ===========================================================
    # Load data
    load('Haim et al JOP Replication - Example Nets.RData')
      # Object: 'low' - an igraph barangay network in our sample with low embeddedness
      # Object: 'high' - an igraph barangay network in our sample with high embeddedness
 
    # layout
    set.seed(1234)
    layout.low<-layout.fruchterman.reingold(low)
    layout.high<-layout.fruchterman.reingold(high)

    # Define Citizen Size
    # Reduce skew in degree so all nodes are visible
    V(low)$degree[V(low)$degree>30]<-V(low)$degree[V(low)$degree>30]-15
    V(low)$size<-V(low)$degree/8+4
    V(high)$degree[V(high)$degree>30]<-V(high)$degree[V(high)$degree>30]-5
    V(high)$size<-V(high)$degree/8+4

    # Define tanod size
    V(low)$size[V(low)$tanod==1]<-(V(low)$size[V(low)$tanod==1]+1.5)*1.3
    V(high)$size[V(high)$tanod==1]<-(V(high)$size[V(high)$tanod==1]+1.5)*1.3    

    # Define shapes
    V(low)$shape<-'circle'
    V(low)$shape[V(low)$tanod==1]<-'square'    
    V(high)$shape<-'circle'
    V(high)$shape[V(high)$tanod==1]<-'square'    

    # Color Nodes by social distance
    V(low)$color<-brewer.pal(9, "Reds")[1]      
    V(low)$color[V(low)$connect.min==1]<-brewer.pal(9, "Reds")[6]
    V(low)$color[V(low)$connect.min==2]<-brewer.pal(9, "Reds")[3]
    V(low)$color[V(low)$tanod==1]<-brewer.pal(9, "Reds")[9]

    V(high)$color<-brewer.pal(9, "Reds")[1]            
    V(high)$color[V(high)$connect.min==1]<-brewer.pal(9, "Reds")[6]
    V(high)$color[V(high)$connect.min==2]<-brewer.pal(9, "Reds")[3]
    V(high)$color[V(high)$tanod==1]<-brewer.pal(9, "Reds")[9]

  # ===========================================================
  # Create Plots
  # ===========================================================
    # Low Embeddedness Plot
      par(mar=c(1,1,1,1),bg=NA)
      plot(low, layout=layout.low, vertex.label=NA, vertex.size=V(low)$size, vertex.shape=V(low)$shape, edge.color=adjustcolor('grey',.85),edge.curved=.2, edge.width=1.4)

    # High Embeddedness Plot
      par(mar=c(1,1,1,1),bg=NA)
      plot(high, layout=layout.high, vertex.label=NA, vertex.size=V(high)$size, vertex.shape=V(high)$shape, edge.color=adjustcolor('grey',.85),edge.curved=.2, edge.width=1.4)

    # NOTE: Legend was added in Preview


#################################################
######### Figure 2: Dyadic Descriptives #########
#################################################

  # Load Dyadic survey data
  dyads<-as.data.table(read.dta13('Haim et al JOP Replication - Dyadic Level.dta'))

  # Create a datatable for the 3 questions
    trust<-c(mean(na.omit(dyads[geodist==1]$tanod_trust)),
      mean(na.omit(dyads[geodist==2]$tanod_trust)),
      mean(na.omit(dyads[geodist==3]$tanod_trust)),
      mean(na.omit(dyads[geodist==4]$tanod_trust)),
      mean(na.omit(dyads[geodist>4]$tanod_trust)))
    report<-c(mean(na.omit(dyads[geodist==1]$tanod_report)),
      mean(na.omit(dyads[geodist==2]$tanod_report)),
      mean(na.omit(dyads[geodist==3]$tanod_report)),
      mean(na.omit(dyads[geodist==4]$tanod_report)),
      mean(na.omit(dyads[geodist>4]$tanod_report)))
    fair<-c(mean(na.omit(dyads[geodist==1]$tanod_fair)),
      mean(na.omit(dyads[geodist==2]$tanod_fair)),
      mean(na.omit(dyads[geodist==3]$tanod_fair)),
      mean(na.omit(dyads[geodist==4]$tanod_fair)),
      mean(na.omit(dyads[geodist>4]$tanod_fair)))    
    summary<-data.table(geodist=c('1','2','3','4','5+'),trust=trust,report=report,fair=fair)

  # --- Create Barplots --- #
    p.trust<-ggplot(data=summary, aes(x=geodist, y=trust)) +
      geom_bar(stat="identity", fill="steelblue")+
      labs(x='Social Distance',y='',title='Trust')+
      coord_cartesian(ylim=c(3.0,3.6))+
      theme_minimal()+
      theme(plot.margin = unit(c(1,.2,.2,.2), "cm"),axis.title=element_text(size=14),legend.text=element_text(size=14),legend.title=element_text(face='bold',size=14),axis.text.x=element_text(size=14),axis.text.y=element_text(size=14),plot.title = element_text(hjust = 0.5,face='bold',size=16)) 
    p.report<-ggplot(data=summary, aes(x=geodist, y=report)) +
      geom_bar(stat="identity", fill="steelblue")+
      labs(x='Social Distance',y='',title='Report')+
      coord_cartesian(ylim=c(3.0,3.6))+
      theme_minimal()+
      theme(plot.margin = unit(c(1,.2,.2,.2), "cm"),axis.title=element_text(size=14),legend.text=element_text(size=14),legend.title=element_text(face='bold',size=14),axis.text.x=element_text(size=14),axis.text.y=element_text(size=14),plot.title = element_text(hjust = 0.5,face='bold',size=16)) 
    p.fair<-ggplot(data=summary, aes(x=geodist, y=fair)) +
      geom_bar(stat="identity", fill="steelblue")+
      labs(x='Social Distance',y='',title='Fair')+
      coord_cartesian(ylim=c(3.0,3.6))+
      theme_minimal()+
      theme(plot.margin = unit(c(1,.2,.2,.2), "cm"),axis.title=element_text(size=14),legend.text=element_text(size=14),legend.title=element_text(face='bold',size=14),axis.text.x=element_text(size=14),axis.text.y=element_text(size=14),plot.title = element_text(hjust = 0.5,face='bold',size=16)) 
    grid.arrange(p.trust, p.report, p.fair, ncol=3)


###########################################
######### Figure 4: Tanod Balance #########
###########################################

  # Load tanod survey
  tsvy<-fread('Haim et al JOP Replication - Tanod Survey.csv')

  # Dichotomize tanod centrality at the median
  tsvy[,degree.hl:=ifelse(degree<median(na.omit(tsvy$degree)),'Low','High')]

  # Generate scatterplots
  age<-ggerrorplot(tsvy[!is.na(degree.hl)], x = c("degree.hl"), y = c('age'), size=1.5, desc_stat = "mean_sd", add = "jitter", add.params = list(color = "tomato"), ylab='Years',xlab='',font.x=c(20,'bold'),font.y=c(20,'bold'),font.xtickslab=c(18),font.ytickslab=c(18),font.title=c('bold',20),ylim=c(0,80),title='Age')  + theme(plot.title = element_text(hjust = 0.5), axis.title.x = element_text(vjust=-.5))
  years.resident<-ggerrorplot(tsvy[!is.na(degree.hl)], x = c("degree.hl"), y = c('years.resident'), size=1.5, desc_stat = "mean_sd", add = "jitter", add.params = list(color = "royalblue1"), ylab='',xlab='Tanod Centrality',font.x=c(20,'bold'),font.y=c(16,'bold'),font.xtickslab=c(18),font.ytickslab=c(18),font.title=c('bold',20),ylim=c(0,80),title='Resided') + theme(plot.title = element_text(hjust = 0.5), axis.title.x = element_text(vjust=-.5), axis.line.y=element_blank(),axis.text.y=element_blank(),
    axis.ticks.y=element_blank())
  years.served<-ggerrorplot(tsvy[!is.na(degree.hl)], x = c("degree.hl"), y = c('years.served'), size=1.5, desc_stat = "mean_sd", add = "jitter", add.params = list(color = "green3"), ylab='',xlab='',font.x=c(20,'bold'),font.y=c(16,'bold'),font.xtickslab=c(18),font.ytickslab=c(18),font.title=c('bold',20),ylim=c(0,80),title='Served') + theme(plot.title = element_text(hjust = 0.5), axis.title.x = element_text(vjust=-.5),axis.line.y=element_blank(),axis.text.y=element_blank(),
    axis.ticks.y=element_blank())

  # Plot
  grid.arrange(age,years.resident,years.served, ncol=3)

















